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Abstract.  Q-ball  imaging  (QBI)  is  a  high  angular  resolution  diffusion  imaging 
(HARDI)  technique  which  has  been  proven  very  successful  in  resolving 
multiple  intravoxel  fiber  orientations  in  MR  images.  The  standard  computation 
of  the  orientation  distribution  function  (ODE,  the  probability  of  diffusion  in  a 
given  direction)  from  q-ball  uses  linear  radial  projection,  neglecting  the  change 
in  the  volume  element  along  the  ray,  thereby  resulting  in  distributions  different 
from  the  true  ODEs.  A  new  technique  has  been  recently  proposed  that,  by 
considering  the  solid  angle  factor,  uses  the  mathematically  correct  definition  of 
the  ODF  and  results  in  a  dimensionless  and  normalized  ODF  expression  from  a 
single  q-shell.  In  this  paper,  we  extend  this  technique  in  order  to  exploit 
HARDI  data  from  multiple  q-shells.  We  consider  the  more  flexible  multi¬ 
exponential  model  for  the  diffusion  signal,  and  show  how  to  efficiently 
compute  the  ODFs  in  constant  solid  angle.  We  describe  our  method  and 
demonstrate  its  improved  performance  on  both  artificial  and  real  HARDI  data. 

Keywords:  Multiple  q-shells,  multi-exponential  models,  orientation 
distribution  function  (ODF),  q-ball  imaging  (QBI),  high  angular  resolution 
diffusion  imaging  (HARDI),  diffusion-weighted  magnetic  resonance  imaging 
(DWMRI),  solid  angle. 


1  Introduction 

Diffusion-weighted  magnetic  resonance  imaging  (DWMRI)  provides  valuable 
information  about  the  fiber  architecture  of  neural  tissue  by  measuring  the  diffusion  of 
water  molecules  in  three-dimensional  (3D)  space.  The  diffusion  function  may  be 
measured  by  using  the  model-free  diffusion  spectrum  imaging  (DSI)  [1],  which  is  the 
direct  Fourier  inversion  of  the  diffusion  signal.  This  technique  is  however  time 
intensive,  as  it  measures  the  diffusion  signal  on  a  3D  Cartesian  lattice.  Thus,  an 
alternative  approach  based  on  sampling  on  one  or  multiple  spherical  shells  has  been 
proposed,  referred  to  as  high  angular  resolution  diffusion  imaging  (HARDI)  [2]. 

While  the  3D  probability  density  function  (PDF)  of  the  diffusion  is  helpful  in 
studying  the  tissue  microstructure,  the  orientation  distribution  function  (ODF)  -  the 
marginal  probability  of  diffusion  in  a  given  direction  -  is  the  quantity  of  interest  for 
mapping  the  orientation  architecture  of  the  tissue.  Q-ball  imaging  (QBI),  [3],  is  a 
widely  used  ODF  reconstruction  scheme  for  HARDI,  based  on  a  spherical 
tomographic  inversion  called  the  Funk-Radon  transform.  This  technique’s  simplicity 


and  its  ability  to  resolve  intravoxel  fiber  orientations  have  made  it  popular  for  fiber 
tracking  and  characterizing  white  matter  architecture.  Moreover,  a  few  works  have 
suggested  exploiting  data  from  multiple  q-shells  to  benefit  from  the  high  signal-to- 
noise  ratio  (SNR)  and  high  angular  contrast-to-noise  ratio  (CNR)  of  the  data  acquired 
at  respectively  low  and  high  b-values,  [3]-[5].  Using  multiple  q-shells  also  allows  us 
to  employ  richer  models  for  the  diffusion  signal,  as  discussed  in  this  paper. 

Nonetheless,  with  the  exception  of  our  previous  paper  [6]  and  a  very  recent  parallel 
and  independent  work  [7]  (the  differences  will  he  detailed  in  Sec.  2.2),  the  definition 
of  the  ODF  used  in  QBI  has  been  different  from  the  actual  marginal  PDF  of  diffusion 
in  a  constant  solid  angle.  It  has  been  computed  as  a  linear  radial  projection  of  the 
PDF,  which  does  not  take  into  account  the  quadratic  growth  of  the  volume  element 
with  respect  to  its  distance  from  the  origin  (see  Sec.  2. 1  for  details).  This  inaccurate 
formulation  generally  distorts  the  ODF,  and  has  created  the  need  for  post-processing 
such  as  manual  normalization  and  sharpening  [8]. 

We  recently  proposed,  [6],  a  new  ODF  expression  for  QBI  which  is  derived  from 
the  proper  definition  of  the  ODF  in  constant  solid  angle.  We  showed  that  the 
computed  ODF  is  inherently  normalized  and  dimensionless,  producing  without  any 
post-processing,  sharp  ODFs  with  improved  resolution  of  multiple  fiber  orientations. 
In  this  paper,  we  extend  this  work  by  deriving  a  general  formulation  for  multiple  q- 
shell  QBI.  We  demonstrate  the  improvement  achieved  by  considering  the  information 
from  multiple  q-shells,  and  using  richer  multi-exponential  models. 

In  Sec.  2  we  describe  the  foundation  of  our  mathematical  derivation,  along  with  a 
brief  version  of  the  proof,  and  also  provide  an  implementation  scheme.  Experimental 
results  are  presented  in  Sec.  3,  along  with  a  brief  discussion. 


2  ODF  Computation  in  Solid  Angle:  Multiple  q-Shell  Formulation 

2.1  General  ODF  Definition 

The  PDF  of  the  diffusion  of  water  molecules,  P(r),  gives  the  displacement 
probability  P{r)dv  of  a  molecule,  initially  placed  at  the  origin,  to  be  in  the 
infinitesimal  volume  dv  located  at  r  after  a  certain  amount  of  time.  We  assume  this 
function  to  be  symmetric  (i.e.  P(— r)  =  P(r)),  which  is  a  quite  common  assumption 
in  DWMRI.  The  PDF  is  represented  in  the  standard  spherical  coordinates,  (r,  6,  0), 
with  the  displacement  vector  r  =  ru  and  the  unit  direction  vector  u(6, 0)  = 
(sin  6  cos  0 ,  sin  S  sin  0 ,  cos  6)^.  The  volume  element  in  this  case  is  dv  =  r^drdQ. 
with  dQ.  =  sin  9  dOdcj)  being  the  infinitesimal  solid  angle  element. 

We  denote  by  ODF{u)dO.  the  probability  of  diffusion  in  the  direction  u  through 
the  solid  angle  dPl,  which  is  computed  by  integrating  the  displacement  probabilities, 
i.e.,  P{f}dv  =  P(ru)r^drdQ.,  for  all  magnitude  r,  while  keeping  u  constant: 

J-  00 

P(ru)r^dr  (1) 

0 

The  above  definition,  which  is  normalized  and  dimensionless,  is  the  integral  of  the 
probability  values  in  a  cone  of  “very  small”  constant  solid  angle.  This  correct 


definition  was  used  for  instance  by  the  authors  of  [1]  in  DSI,  where  P(r)  was  first 
computed  from  the  diffusion  data  via  Fourier  inversion  and  then  integrated  to 
calculate  the  ODF.  However  to  the  best  of  our  knowledge,  the  expression  for  ODF 
reconstruction  so  far  used  in  QBI  [3],  (except  for  our  previous  work  [6]  and  a  very 
recent  parallel  and  independent  paper  [7],  both  for  single  q-shell)  is  different  from  Eq. 
(1),  in  the  sense  that  the  integral  is  not  weighted  by  the  important  (and  mathematically 
correct)  factor  Without  including  this  factor,  the  radial  projection  gives  an 
artificial  weight  to  P(r)  which  is,  respectively,  too  large  and  too  small  for  points 
close  to  and  far  from  the  origin.  Moreover,  the  ODF  will  not  be  necessarily 
normalized  or  dimensionless,  and  manual  normalization  will  be  required. 

Next  we  derive  a  closed-form  ODF  expression  in  multiple  q-shell  QBI  using  the 
correct  -weighted  integral. 


2.2  Q-ball  Imaging  ODF  Reconstruction 

In  this  section,  we  derive  the  ODF  expression  in  multiple  q-shell  QBI,  and  present  a 
brief  proof  of  the  derivation. 

Let  E{fl)  be  the  3D  Fourier  transform  function  of  P{f).  Theoretically,  we  know 
that  ECO)  =  1,  since  the  zero  frequency  of  a  PDF  is  its  integral  over  the  entire  space, 
yielding  1.  In  addition,  we  have  the  values  of  E{q)  measured  on  M  different  q-balls, 
i.e.,  the  frequencies  with  constant  norm  |q|  =  q^,  i  =  1,  as  E^iu)  ■=  E^q^u)  = 

where  Si(u)  is  the  HARDI  signal  on  the  t*  q-ball  and  is  the  base-line  image. 

Our  mathematical  derivation  is  based  on  the  following  two  relatively  simple  yet 
fundamental  facts  from  Fourier  analysis: 

•  The  Fourier  transform  of  P(r)|r|^  is  —V^E(jq),  where  is  the  Laplacian 
operator. 

•  For  a  symmetric  function  /:  ^  M  with  the  3D  Fourier  transform  function 

F(q),  and  for  the  arbitrary  unit  vector  u,  we  have  that  f  (ru)dr  = 

with  u-'-  being  the  plane  perpendicular  to  u. 

Combining  these  statements  with  Eq.  (1)  leads  to 

ODFiu)  =  -  ^  II^^V^E(q)d^q 

Now,  without  loss  of  generality,  we  choose  our  coordinates  such  that  z  =  u,  thus 
making  u-‘-  the  qx-qy  plane.  We  then  use  the  following  expansion  for  the  Laplacian  in 
spherical  coordinates,  (q,  6, 0): 

,  1  ^  1  , 

^  =  n 

q  aq'^  q'^ 

where  is  the  Laplace-Beltrami  operator,  which  is  defined  independently  of  the 
radial  component  q,  as  ViE  =  — -^fsinP— 1  -I - The  surface  integral  on 

r  't’  D  Sinsaev  36 J  sin^  e  302  & 

the  qx-qy  plane  is  computed  by  fixing  9  =  -  and  using  the  expression  d^q  =  qdqdcj) 
for  the  surface  element 


ODF(z)  =  - 


r^n  r'i 

5^1  i 


^^E{q)qdqd(p 

32 


qdq  d<p 


The  integral  of  the  first  term  can  be  seen  to  be  constant  and  independent  of  E{q), 

r2n  r“  /p  g2  \ 

—  (qE)  jqdqdcj)  =  —2nE(0)  =  —2n 


rf 

*'0  “'0 


q  dq^ 


Therefore, 


/ 


n0( 

) 


VlE(q)dq  d(p 


while  6  =  -  is  kept  constant  in  the  integration. 

To  compute  the  integral  of  the  second  term,  the  values  of  E{q)  are  required  in  the 
entire  q-space,  which  are  in  general  -  except  for  the  time-consuming  DSI  modality  - 
not  available.  Thus,  we  need  to  approximate  E(q)  from  the  values  measured  on  the  q- 
balls.  In  this  work,  we  consider  the  following  radial  multi-exponential  model  [9], 


E(qu) 


with  the  constraints 


0  <  afc(ti),Afc(u)  <  1 
N 

^4(u)  =  l  (2) 

k=l 

where  Eq.  (2)  comes  from  the  fact  that  E{0)  =  1.'  Once  the  values  of  2,^  and  are 
estimated  (see  Sec.  2.3),  they  can  be  used  in  the  following  ODE  expression,  which  is 
obtained  by  a  few  more  steps  of  calculation, 

1  1 

ODF(z)  =  ^  +  J  Vg  Ifc(u)  ln(-  In  a^(u))  d(t) 

°  fe=i 

Finally,  rewriting  the  expression  independent  of  the  choice  of  the  axes,  the 
following  analytical  formula  can  be  derived  for  the  ODE: 

{N 

J^Afc(u)ln(-lnafc(u))  (3) 

where  Fl?r(/(u)}  ;= /(w)5(|w|  —  l)d^w  is  the  Funk-Radon  transform  [3], 
with  5(0  being  the  Dirac  delta  function. 

The  above  ODE  expression  is  dimensionless  and  intrinsically  normalized,  since  the 
integrals  of  the  first  and  second  terms  over  the  sphere  are  respectively  1  and  0.  This  is 
in  contrast  to  the  (single  q-shell)  ODE  formulas  used  in  original  QBI,  i.e., 
^FRT{Ej^(u)],  and  also  in  [7],  where  a  normalization  factor  Z  is  needed.  Additional 
differences  can  be  observed  in  the  approach  presented  here  and  in  [6],  compared  to 
[7].  As  demonstrated  here,  integration  of  the  radial  part  of  the  Laplacian  on  the  plane 


'  This  is  in  fact  an  additional  advantage  of  this  model  over  the  original  QBI  model  for  single  q- 
shell,  i.e.  Eiqu)  =  E{qiU)S{q  —  q^),  where  £(0)  was  assumed  to  be  zero. 


always  results  in  a  constant  without  requiring  any  model  for  the  diffusion  signal.  Yet, 
[7]  uses  the  Bessel  approximation  of  the  Dirac  delta  function  which  yields  a  variable 
(sometimes  negative)  term.  As  for  the  integral  of  the  tangential  term  of  the  Laplacian, 
we  use  the  exponential  model  that  is  particularly  consistent  with  ^CO)  =  1,  in  contrast 
to  [7]  that  assumes  the  tangential  term  to  be  zero  outside  the  q-ball,  leading  to  an 
expression  similar  to  Laplacian-Beltrami  sharpening. 


2.3  Parameter  Estimation 

In  order  to  approximate  the  diffusion  signal  in  a  direction  u  by  a  weighted  sum  of  N 
exponentials,  we  need  to  estimate  the  2N  parameters  Afc(u)  and  for  k  = 

1, ... ,  N.  We  continue  this  subsection  considering  a  fixed  direction,  and  therefore  drop 
the  notation  (u).  To  estimate  the  aforementioned  parameters,  at  least  2N  —  1 
independent  equations  -  besides  Eq.  (2)  -  are  required,  which  can  be  obtained  from 
the  HARDI  signals  measured  on  M  q-balls,  for  M  >  2N  —  1.  Numerical  optimization 
approaches  such  as  the  trust  region  algorithm,  [10],  may  be  employed  to  solve  this 
non-linear  system  in  the  most  general  case.  Here,  however,  we  discuss  two  special 
cases  with  closed-form  analytical  solutions. 

The  mono-exponential  assumption  (N  =  1)  requires  measurement  on  at  least 

1 

M  =  1  q-ball.  As  it  has  been  shown  in  [6],  M  =  1  leads  to  2,]^  =  1  and 
Furthermore,  if  measured  values  are  provided  on  more  than  one  q-balls  and  the  mono¬ 
exponential  model  is  still  desired,  one  can  fit  the  best  exponential  by  computing  the 

1  s 

average  Apparent  Diffusion  Coefficient  {ADC  ■=  — In—)  across  all  the  q-shells. 

b  So 

Another  practical  case  of  great  interest  arises  when  we  consider  the  richer  bi¬ 
exponential  model  (iV  =  2,  see  for  example  [11])  to  reconstruct  the  ODFs  from  (at 
least)  M  =  3  q-shells.  Parameterizing  the  problem  in  terms  of  b-values,  h;  =  rqf,  and 
choosing  the  physical  units  such  that  the  diffusion  time  t  =  1  (see  also  Footnote  2), 
we  obtain  (for  M  =  3)  the  following  system  of  equations  for  each  direction: 

+  (1  -  X)P’’i  =Ei  ,  i  =  1,2,3 
0  <  a,p,X  <  1 

An  analytical  solution  can  be  derived  for  the  particular  and  reasonable  case  when 
the  sequence  0,bj^,b2,b^  is  an  arithmetic  progress.^  We  describe  this  solution  here, 
along  with  some  regularization  that  guarantees  the  parameters  to  remain  within  the 
correct  range."*  Without  loss  of  generality,  let  us  assume  a  >  /?,  and  also  choose  the 
physical  units  such  that  b^  =  1,  hj  =  2,  and  b^  =  3.  Then, 

la*  +  (1  -  2)/?‘  =  Ei  ,  i  =  1,2,3 

We  first  define  and  calculate  the  following  two  quantities: 


^  Note  that  if  the  set  {af.{u)}  is  a  solution,  then  {^(u)’'}  for  a  constant  y  can  be  shown  to 
result  in  the  same  computed  ODF.  Therefore,  since  in  the  mono-exponential  case  y  ;=  qf  is  a 
constant,  (u)  =  (u)  is  also  a  correct  solution. 

^  The  sequence  X2, ... ,  Xj, ...  is  an  arithmetic  progress  if  x,  —  Xi_i  is  constant. 

"*  Recall  that  the  three  parameters  can  also  be  computed  in  the  general  case  following 
optimization  techniques  such  as  those  in  [9].  The  proposed  ODF  model  is  general,  and  it 
becomes  only  simpler  when  the  data  is  acquired  at  an  arithmetic  sequence  of  b-values. 


_  a+P  ^  E3  _  a _ ^  ^  /  E-i  £1^2  \  _  ^1^3  El 

2  2{E^-El)  ’  '  2  ^[2(E2-Eiy  E^  -  E^ 

The  parameters  are  afterward  computed  as  follows: 

a=A+B  ,  p  =  A-B  ,  ^  =  l  + 

However,  we  still  need  to  ensure  that  they  are  real  and  in  the  correct  ranges.  One 
can  verify  that  these  conditions  are  satisfied  by  enforcing  the  following  constraints: 

0  <  E3  <  ^2  <  ^1  <  1  ,  Ef  <  E2  I  E2  <  E-^E^ 

£*3  —  E1E2  ^  E2  £"1  +  E^E^  ~  £2 

Thus,  we  can  obtain  the  optimal  values  of  a,  p,  and  1,  by  initially  projecting  £iS 
onto  the  subspace  defined  by  these  inequalities,^  and  then  computing  the  parameters. 

2.4  Implementation 

Our  implementation  of  the  ODF  reconstruction  from  the  estimated  values  of  Afe(ti) 
and  afc(ti)  makes  use  of  the  spherical  harmonic  (SH)  basis,  Yl^{u),  which  is  common 
for  the  analysis  of  HARDI  data.  The  steps  taken  here  to  numerically  compute  Eq.  (3) 
are  similar  to  those  described  in  [8].  Particularly,  we  use  the  real  and  symmetric 
modified  SH  basis  introduced  in  [8],  where  SH  functions  are  indexed  by  a  single 
parameter  j  corresponding  to  Ij  and  m,-.  We  adopt  a  minimum  least  square  scheme  to 
compute  a  set  of  modified  SH  coefficients,  Cj,  such  that 

N  R 

^4(w)ln(-lnafe(tt))  «  I  CjYj(u) 

k=l  j=l 

where  R  =  {L  +  1)(L  +  2)/2  ,  with  L  being  the  order  of  the  SH  basis  (we  chose 
L=  4  throughout  our  experiments).  Next,  since  the  SH  elements  are  eigenfunctions  of 
the  Laplace-Beltrami  operator,  we  compute  VlY.k=i2kiE.)\n(—\nai^(u))  by 
multiplying  the  coefficients  Cj  by  their  corresponding  eigenvalues,  -  lj{lj  +  l).  Then, 
as  suggested  in  [8],  the  Funk-Radon  transform  is  computed  by  multiplying  the 
coefficients  by  2nPi.(fY),  where  £;(•)  is  the  Legendre  polynomial  of  degree  I,  with 

^i(O)  =  fot  even  1.  Finally,  given  that  T^(u)=— the  SH 

2x4'X'**xf  '2.'\Jtc 

coefficients  of  the  ODF  are  derived  as 


r^^“"^^2x4x...x0,-2)^^ 

By  taking  advantage  of  the  SH  framework,  the  implementation  of  the  proposed 
formula  for  the  true  ODF  is  as  straightforward  as  the  one  introduced  in  [8]  for  the 
original  ODF  formula,  or  even  simpler  if  no  further  sharpening  is  to  be  performed. 


1 


X  3  X  •••  X  {ij  +  1) 


^  Note  that  such  projection  is  usually  necessary,  because  the  bi-exponential  assumption  may  not 
be  accurate  and  the  data  may  be  noisy.  Moreover,  using  a  small  separating  margin  of 
S  =  0.01~0.1  in  the  inequalities  makes  the  ODFs  in  practice  more  stable. 


3  Results  and  Discussion 


To  demonstrate  the  advantages  of  exploiting  multiple  q-shells  in  QBI,  we  first  show 
the  experimental  results  on  an  artificial  example  which  consists  of  large  diffusion 
values  in  two  orthogonal  directions.  We  synthesized  diffusion  images  by  sampling  the 

sum  of  two  exponentials,  E^q)  =  +  ^|cos  (p]’  )  seven  q-shells 

(b  =  =  1,2,  ...,7)  and  in  76  directions,  uniformly  distributed  on  the  hemisphere. 

Figure  1  illustrates  the  ODFs  reconstructed  from  single  q-shells  for  different  b-values, 
three  q-shells  with  mono-exponential  model,  and  three  q-shells  with  bi-exponential 
model.  As  can  be  observed,  for  the  data  acquired  at  low  b-values  (b  =  1,2,3),  the  bi¬ 
exponential  model  using  three  q-shells  is  the  only  method  correctly  resolving  the 
horizontal  and  vertical  ODF  peaks,  corresponding  to  the  strong  ADC  values  in  those 
directions  (0  =  0°,  90°,  180°,  270°).  It  should  be  noted,  however,  that  the  drawback 
of  such  a  more  general  model  is  its  lesser  robustness  to  noise,  as  low  order  models  are 
often  more  robust  (e.g.,  computing  the  average  of  a  signal  is  more  robust  than 
estimating  the  actual  signal).  ODFs  are  shown  as  they  are;  no  min-max  normalization 
is  used  in  any  of  the  figures.  Dark  red  represents  negative  values. 

We  also  tested  our  method  on  the  real  HARDI  dataset  introduced  in  [12].  An 
anesthetized  young  Macaca  mulatta  monkey  was  scanned  using  a  7T  MR  scanner 
(Siemens)  equipped  with  a  head  gradient  coil  (80mT/m  G-maximum,  200mT/m/ms) 
with  a  diffusion  weighted  spin-echo  EPI  sequence.  Diffusion  images  were  acquired 
(twice  during  the  same  session,  and  then  averaged)  over  100  directions  uniformly 
distributed  on  the  sphere.  We  used  three  b-values  of  1000,  2000,  and  3000  s/mm^, 
TR/TE  of  4600/65  ms,  and  the  voxel  size  of  1x1x1  mm^.  The  ODEs  were 
reconstructed  from  the  three  q-shells  using  both  mono-exponential  and  bi-exponential 
methods,  and  also  from  the  single  q-shells  individually.  Figure  2  depicts  the  results  on 
a  coronal  slice  through  the  centrum  semiovale  area,  superimposed  on  the  fractional 
anisotropy  (FA)  map.  Note  how  using  the  bi-exponential  method  allows  for  more 
clear  recovery  of  certain  fiber  bundles,  such  as  callosal  radiations  and  corticospinal 
tract,  and  better  resolution  of  crossing  areas  (see  outlined  regions  in  Fig.  2).  Figure  2 
(top,  right)  is  the  only  subfigure  illustrating  results  by  the  original  QBI  (without  r^). 
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